#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <float.h>
#include <limits.h>
#include <assert.h>
#include <mpi.h>

#include "halofinder.h"

#define mabs(A) ((A) < 0.0 ? -(A) : (A))

typedef struct {
  float aexpn;
  float aexp0;
  float amplt;
  float astep;
  int   istep;
  float partw;
  float tintg;
  float ekin;
  float ekin1;
  float ekin2;
  float au0;
  float aeu0;
  int   Nrow;
  int   Ngrid;
  int   Nspecies;
  int   Nseed;
  float Om0;
  float Oml0;
  float hubble;
  float Wp5;
  float Ocurv;
  float Omb0;
  float mass[10];
  int   num[10];
  float fill[80];
} particle_header;

typedef struct {
  float aexpn;
  float aexp0;
  float amplt;
  float astep;
  int   istep;
  float partw;
  float tintg;
  float ekin;
  float ekin1;
  float ekin2;
  float au0;
  float aeu0;
  int   Nrow;
  int   Ngrid;
  int   Nspecies;
  int   Nseed;
  float Om0;
  float Oml0;
  float hubble;
  float Wp5;
  float Ocurv;
  float mass[10];
  int   num[10];
  float fill[80];
} nbody_particle_header;

void reorder( char *buffer, size_t size ) {
  size_t i;
  char tmp;

  for ( i = 0; i < (size/2); i++ ) {
    tmp = buffer[i];
    buffer[i] = buffer[size - i - 1];
    buffer[size - i - 1] = tmp;
  }
}

int read_ART_Page(long PageNum, long *NumBuffParts, float **buffPx, float **buffPy, float **buffPz, float **buffVx, float **buffVy, float **buffVz)
{
  int i,j;//, nCnt;
  nbody_particle_header nbody_header;
  
  FILE *header_file;
  FILE *data_file;
  int size;
  int nbody_flag;
  int endian;
  int num_grid;
  int num_particles;
  int page_size;
  int num_pages;	
  //int num[10];
  //float mass[10];
  //float *input_page;
  int num_parts_in_page;
  char desc[45];
  //long baseInd;
  
  float vg2kms;
  float box_size;
  //int jcnt;//,rcube;//,rhalf;
  //float xx[6];
  
  
  char header[1024];
  char data[1024];
  
  header_file = fopen(haloFinderData.SnapshotHeaderList,"r");
  assert(header_file != NULL);
  fscanf(header_file,"%s\n",header);
  fclose(header_file);
  
  data_file = fopen(haloFinderData.SnapshotFileList,"r");
  assert(data_file != NULL);
  fscanf(data_file,"%s\n",data);
  fclose(data_file);
  
  box_size = (float) (haloFinderData.BoxLengthSnapshot);
  
  nbody_flag = 0;
  endian = 0;
  
  header_file = fopen(header, "r");
  if ( header_file == NULL ) {
    fprintf(stderr,"Unable to open %s!\n", header );
    exit(1);
  }
  
  data_file = fopen(data,"rb");
  if ( data_file == NULL ) {
    fprintf(stderr,"Unable to open %s!\n", data );
    MPI_Abort(MPI_COMM_WORLD,123);
  }
  
  fread( &size, sizeof(int), (size_t) 1, header_file );
  
  if ( size == sizeof(nbody_particle_header)+45 ) {
    nbody_flag = 1;
  } else if ( size != sizeof(particle_header)+45 ) {
    reorder( (char *)&size, sizeof(int) );
        
    if ( size == sizeof(particle_header)+45 ) {
      endian = 1;
    } else if ( size == sizeof(nbody_particle_header)+45 ) {
      endian = 1;
      nbody_flag = 1;
    } else {
      fprintf(stderr,"Error reading from particle header, size = %d!\n", size );
      MPI_Abort(MPI_COMM_WORLD,123);
    }
  }
  
  if( nbody_flag == 1)
    nbody_flag = 1;
  
  fread( desc, sizeof(char), (size_t) 45, header_file );
  
  fread( &nbody_header, sizeof(nbody_particle_header), (size_t) 1, header_file );
  
  if ( endian ) {
    reorder( (char *)&nbody_header.Ngrid, sizeof(int) );
    reorder( (char *)&nbody_header.Nrow, sizeof(int) );
    reorder( (char *)&nbody_header.Nspecies, sizeof(int) );
    reorder( (char *)&nbody_header.aexpn, sizeof(float) );
    reorder( (char *)&nbody_header.Om0, sizeof(float) );
    reorder( (char *)&nbody_header.Oml0, sizeof(float) );
    reorder( (char *)&nbody_header.hubble, sizeof(float) );
    reorder( (char *)&nbody_header.istep, sizeof(int) );
    
    for ( i = 0; i < nbody_header.Nspecies; i++ ) {
      reorder( (char *)&nbody_header.num[i], sizeof(int) );
      reorder( (char *)&nbody_header.mass[i], sizeof(float) );
    }
  }
  
  /*for ( i = 0; i < nbody_header.Nspecies; i++ ) {
    num[i] = nbody_header.num[i];
    mass[i] = nbody_header.mass[i];
    }
  */
  
  num_particles = nbody_header.num[ nbody_header.Nspecies - 1];
  
  //fprintf(stderr,"Nrow = %u\n", nbody_header.Nrow );
  page_size = nbody_header.Nrow*nbody_header.Nrow;
  num_grid = nbody_header.Ngrid;
  
  
  fclose( header_file );
  
  num_pages = (num_particles-1) / page_size + 1;
  
  // fprintf(stderr,"endian = %u\n", endian );
  // fprintf(stderr,"num_grid = %u\n", num_grid );
  //fprintf(stderr,"%d particles\n", num_particles );
  // fprintf(stderr,"Page size = %u\n", page_size );
  // fprintf(stderr,"%u pages\n", num_pages );
  
  //rcube = box_size;
  //rhalf = rcube/2;
  
  /* allocate page */
  //input_page = malloc( 6*page_size * sizeof(float) );
  //assert(input_page != NULL);
	 
  //*np = num_particles*particle_fraction;
  
  //fprintf(stderr,"Allocation space for %d partciles (%.2f) %f\n",*np,((float) *np)*1./((float) num_particles),particle_fraction);
  
  //*x=(float *)calloc(*np,sizeof(float));
  //*y=(float *)calloc(*np,sizeof(float));
  //*z=(float *)calloc(*np,sizeof(float));
  
  //*vx=(float *)calloc(*np,sizeof(float));
  //*vy=(float *)calloc(*np,sizeof(float));
  //*vz=(float *)calloc(*np,sizeof(float));
 
  //if(*x == NULL || *y == NULL || *z == NULL || *vx == NULL || *vy == NULL || *vz == NULL)
  //{
  //  fprintf(stderr,"out of memory at *x/*y/*z/*vx/*vy/*vz in read_ART...\n");
  //  MPI_Abort(MPI_COMM_WORLD,123);
  //}
 
  //jcnt = 0;
  
  //      vg2kms = 100.0*boxh0 / (aexpn*ngridc)
  vg2kms = 100*box_size/num_grid;
  // fprintf(stderr,"%e vg2kms\n", vg2kms);
  //fprintf(stderr,"z = %e, istep = %d\n",1/nbody_header.aexpn-1,nbody_header.istep);
  //fprintf(stderr,"Omega_m = %f, Omega_lambda = %f, hubble = %f\n",nbody_header.Om0,nbody_header.Oml0,nbody_header.hubble);	 
    
  /*
  ** Read Stuff!
  */
  //nCnt = 0;
  //for (i=0;i<num_pages;++i) {
  i = (int) PageNum;

#ifdef DEBUG  
  if(DEBUG_LEVEL > 1)
    fprintf(stderr,"Page Number %d of %d\n",i,num_pages);
#endif
  if ( i == num_pages - 1 ) {
    num_parts_in_page = num_particles - page_size*(num_pages-1);
  } else {
    num_parts_in_page = page_size;
  }
  *NumBuffParts = (long) num_parts_in_page;
  
  *buffPx = (float*)malloc(sizeof(float)*num_parts_in_page);
  assert((*buffPx) != NULL);
  *buffPy = (float*)malloc(sizeof(float)*num_parts_in_page);
  assert((*buffPy) != NULL);
  *buffPz = (float*)malloc(sizeof(float)*num_parts_in_page);
  assert((*buffPz) != NULL);
  *buffVx = (float*)malloc(sizeof(float)*num_parts_in_page);
  assert((*buffVx) != NULL);
  *buffVy = (float*)malloc(sizeof(float)*num_parts_in_page);
  assert((*buffVy) != NULL);
  *buffVz = (float*)malloc(sizeof(float)*num_parts_in_page);
  assert((*buffVz) != NULL);
  
  //baseInd = PageNum*page_size;
  
  /* easiest to just read in entire page even though we're only using positions */
  if(fseek(data_file,(long) (PageNum*page_size*6*sizeof(float)),SEEK_SET))
    {
      fprintf(stderr,"Error seeking to page %d from data file!\n", i );
      
      MPI_Abort(MPI_COMM_WORLD,123);
      
    }
#ifdef DEBUG  
  if(DEBUG_LEVEL > 1)
    fprintf(stderr,"did seek to page\n");
#endif

  if ( fread(*buffPx, sizeof(float), (size_t) page_size, data_file ) != ((unsigned int) page_size) ) {
    fprintf(stderr,"Error reading page %d from data file! 1\n", i );
    
    MPI_Abort(MPI_COMM_WORLD,123);
  }
  if ( fread(*buffPy, sizeof(float), (size_t) page_size, data_file ) != ((unsigned int) page_size) ) {
    fprintf(stderr,"Error reading page %d from data file! 2\n", i );
    
    MPI_Abort(MPI_COMM_WORLD,123);
  }
  if ( fread(*buffPz, sizeof(float), (size_t) page_size, data_file ) != ((unsigned int) page_size) ) {
    fprintf(stderr,"Error reading page %d from data file! 3\n", i );
    
    MPI_Abort(MPI_COMM_WORLD,123);
  }
  if ( fread(*buffVx, sizeof(float), (size_t) page_size, data_file ) != ((unsigned int) page_size) ) {
    fprintf(stderr,"Error reading page %d from data file! 4\n", i );
    
    MPI_Abort(MPI_COMM_WORLD,123);
  }
  if ( fread(*buffVy, sizeof(float), (size_t) page_size, data_file ) != ((unsigned int) page_size) ) {
    fprintf(stderr,"Error reading page %d from data file! 5\n", i );
    
    MPI_Abort(MPI_COMM_WORLD,123);
  }
  if ( fread(*buffVz, sizeof(float), (size_t) page_size, data_file ) != ((unsigned int) page_size) ) {
    fprintf(stderr,"Error reading page %d from data file! 6\n", i );
    
    MPI_Abort(MPI_COMM_WORLD,123);
  }
#ifdef DEBUG  
  if(DEBUG_LEVEL > 1)
    fprintf(stderr,"read page from disk\n");
#endif
  
  if ( endian ) {
    for ( j = 0; j < num_parts_in_page; j++ ) {
      reorder( (char *)&((*buffPx)[j]), sizeof(float) );
      reorder( (char *)&((*buffPy)[j]), sizeof(float) );
      reorder( (char *)&((*buffPz)[j]), sizeof(float) );
      reorder( (char *)&((*buffVx)[j]), sizeof(float) );
      reorder( (char *)&((*buffVy)[j]), sizeof(float) );
      reorder( (char *)&((*buffVz)[j]), sizeof(float) );
      /*
	reorder( (char *)&input_page[j], sizeof(float) );
	reorder( (char *)&input_page[j+page_size], sizeof(float) );
	reorder( (char *)&input_page[j+2*page_size], sizeof(float) );
	reorder( (char *)&input_page[j+3*page_size], sizeof(float) );
	reorder( (char *)&input_page[j+4*page_size], sizeof(float) );
	reorder( (char *)&input_page[j+5*page_size], sizeof(float) );
      */
    }
  }
#ifdef DEBUG  
  if(DEBUG_LEVEL > 1)
    fprintf(stderr,"reordered page if needed\n");
#endif

  for ( j = 0; j < num_parts_in_page; j++ ) {
    
    (*buffPx)[j] = ((*buffPx)[j]-1)*box_size/num_grid;
    (*buffPy)[j] = ((*buffPy)[j]-1)*box_size/num_grid;
    (*buffPz)[j] = ((*buffPz)[j]-1)*box_size/num_grid;
    (*buffVx)[j] = (*buffVx)[j]*vg2kms;
    (*buffVy)[j] = (*buffVy)[j]*vg2kms;
    (*buffVz)[j] = (*buffVz)[j]*vg2kms;
    //(*BuffParts)[j].index = baseInd + j;
    
    //xx[0] = (input_page[j]-1)*box_size/num_grid;
    //xx[1] = (input_page[j+page_size]-1)*box_size/num_grid;
    //xx[2] = (input_page[j+2*page_size]-1)*box_size/num_grid;
    //xx[3] = input_page[j+3*page_size]*vg2kms;
    //xx[4] = input_page[j+4*page_size]*vg2kms;
    //xx[5] = input_page[j+5*page_size]*vg2kms;
    
    
    //if(drand48()<particle_fraction)
    //{
    //(*x)[jcnt] = xx[0];
    //  (*y)[jcnt] = xx[1];
    //  (*z)[jcnt] = xx[2];
    //  (*vx)[jcnt] = xx[3];
    //  (*vy)[jcnt] = xx[4];
    //  (*vz)[jcnt] = xx[5];
    //  jcnt++;
    //  if(jcnt==*np)goto DONE_READING;
    //}
    //++nCnt;
    
  }
  //}
  /*
    current_specie = 0;
    for ( j = 0; j < num_particles; j++ ) {
    if ( j == num[current_specie] ) {
    current_specie++;
    }
    
    kd->p[j].fMass = mass[current_specie];
    }
  */
  
  //DONE_READING:
  //*np = jcnt;
#ifdef DEBUG  
  if(DEBUG_LEVEL > 1)
    fprintf(stderr,"put parts in buffer\n");
#endif
  
  //free( input_page );
  fclose(data_file);
#ifdef DEBUG  
  if(DEBUG_LEVEL > 1)
    fprintf(stderr,"freed page and closed data file\n");
#endif
  
  return 1;

}

long read_ART_NumPages(void)
{
  int i;
  nbody_particle_header nbody_header;

  FILE *header_file;
  int size;
  int nbody_flag;
  int endian;
  //int num_grid;
  int num_particles;
  int page_size;
  int num_pages;
  char desc[45];
  
  char header[1024];
  header_file = fopen(haloFinderData.SnapshotHeaderList,"r");
  assert(header_file != NULL);
  fscanf(header_file,"%s\n",header);
  fclose(header_file);
  
  nbody_flag = 0;
  endian = 0;

  header_file = fopen(header, "r");
  if ( header_file == NULL )
    {
      fprintf(stderr,"Unable to open %s!\n", header );
      MPI_Abort(MPI_COMM_WORLD,123);
    }

  fread( &size, sizeof(int), (size_t) 1, header_file );

  if ( size == sizeof(nbody_particle_header)+45 )
    nbody_flag = 1;
  else if ( size != sizeof(particle_header)+45 )
    {
      reorder( (char *)&size, sizeof(int) );

      if ( size == sizeof(particle_header)+45 )
        endian = 1;
      else if ( size == sizeof(nbody_particle_header)+45 )
        {
          endian = 1;
          nbody_flag = 1;
        }
      else
        {
          fprintf(stderr,"Error reading from particle header, size = %d!\n", size );
          MPI_Abort(MPI_COMM_WORLD,123);
        }
    }
  
  if( nbody_flag == 1)
    nbody_flag = 1;
  
  fread( desc, sizeof(char), (size_t) 45, header_file );
  fread( &nbody_header, sizeof(nbody_particle_header), (size_t) 1, header_file );
  fclose( header_file );

  if ( endian )
    {
      reorder( (char *)&nbody_header.Ngrid, sizeof(int) );
      reorder( (char *)&nbody_header.Nrow, sizeof(int) );
      reorder( (char *)&nbody_header.Nspecies, sizeof(int) );
      reorder( (char *)&nbody_header.aexpn, sizeof(float) );
      reorder( (char *)&nbody_header.Om0, sizeof(float) );

      for ( i = 0; i < nbody_header.Nspecies; i++ )
        {
          reorder( (char *)&nbody_header.num[i], sizeof(int) );
          reorder( (char *)&nbody_header.mass[i], sizeof(float) );
        }

      reorder( (char *)&nbody_header.aexp0, sizeof(float) );
      reorder( (char *)&nbody_header.amplt, sizeof(float) );
      reorder( (char *)&nbody_header.astep, sizeof(float) );
      reorder( (char *)&nbody_header.istep, sizeof(int) );
      reorder( (char *)&nbody_header.partw, sizeof(float) );
      reorder( (char *)&nbody_header.tintg, sizeof(float) );
      reorder( (char *)&nbody_header.ekin, sizeof(float) );
      reorder( (char *)&nbody_header.ekin1, sizeof(float) );
      reorder( (char *)&nbody_header.ekin2, sizeof(float) );
      reorder( (char *)&nbody_header.au0, sizeof(float) );
      reorder( (char *)&nbody_header.aeu0, sizeof(float) );
      reorder( (char *)&nbody_header.Nseed, sizeof(int) );
      reorder( (char *)&nbody_header.Oml0, sizeof(float) );
      reorder( (char *)&nbody_header.hubble, sizeof(float) );
      reorder( (char *)&nbody_header.Wp5, sizeof(float) );
      reorder( (char *)&nbody_header.Ocurv, sizeof(float) );
      for(i=0;i<80;++i)
        reorder( (char *)&nbody_header.fill[i], sizeof(float) );

    }

  num_particles = nbody_header.num[ nbody_header.Nspecies - 1];
  page_size = nbody_header.Nrow*nbody_header.Nrow;
  //num_grid = nbody_header.Ngrid;
  num_pages = (num_particles-1) / page_size + 1;

  //rcube = box_size;                                                                                                                                                                                   
  //rhalf = rcube/2;                                                                                                                                                                                    
  //vg2kms = 100.0*boxh0 / (aexpn*ngridc)                                                                                                                                                               
  //vg2kms = 100*box_size/num_grid;
  //fprintf(stderr,"Nrow = %u\n", nbody_header.Nrow );
  //fprintf(stderr,"endian = %u\n", endian );
  //fprintf(stderr,"num_grid = %u\n", num_grid );  
  //fprintf(stderr,"Number of particles = %u\n", num_particles );  
  //fprintf(stderr,"Page size = %u\n", page_size );  
  //fprintf(stderr,"Number of pages = %u\n", num_pages );  
  //fprintf(stderr,"vg2kms = %e\n", vg2kms);  
  //fprintf(stderr,"z = %e, istep = %u\n",1/nbody_header.aexpn-1,nbody_header.istep);  
  //fprintf(stderr,"Omega_m = %f\n",nbody_header.Om0);                                                                                                                                                  
  return (long) num_pages;
}

#undef mabs
